Traveling length and minimal traveling time for flow 
through percolation networks with long-range spatial correlations 
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We study the distributions of traveling length I and minimal traveling time t m in through two- 
dimensional percolation porous media characterized by long-range spatial correlations. We model 
the dynamics of fluid displacement by the convective movement of tracer particles driven by a 
pressure difference between two fixed sites ("wells") separated by Euclidean distance r. For strongly 
correlated pore networks at criticality, we find that the probability distribution functions P(l) and 
P{t m in) follow the same scaling Ansatz originally proposed for the uncorrelated case, but with 
quite different scaling exponents. We relate these changes in dynamical behavior to the main 
morphological difference between correlated and uncorrelated clusters, namely, the compactness 
of their backbones. Our simulations reveal that the dynamical scaling exponents di and d t for 
correlated geometries take values intermediate between the uncorrelated and homogeneous limiting 



cases, where I* 
respectively. 



r" 1 and tl 



r and I* and t* min are the most probable values of I and t„ 



PACS numbers: 47.55.Mh, 64.60.Ak 
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I. INTRODUCTION 

Fluid transport in porous media is of central impor- 
tance to problems in petroleum exploration and pro- 
duction Q ||, H Q). The geometry of an oil field can 
be very complex, displaying heterogeneities over a wide 
range of length scales from centimeters to kilometers [|| . 
The most common method of oil recovery is by displace- 
ment. Typically water or a miscible gas (carbon dioxide 
or methane) is injected in a well (or wells) to displace oil 
to other wells. Ultimately the displacing fluid will break 
through into a production well where it must be sepa- 
rated from the oil. At this point, the rate of oil produc- 
tion decreases. For economic purposes, it is important to 
predict when the injected fluid will break through. 

When modeling the process of oil recovery, an open 
question is the effect of the connectedness of the porous 
medium on the dynamical process of fluid displacement. 
If the pore space is so poorly connected as to be con- 
sidered macroscopically heterogeneous, one expects the 
overall behavior of the flowing system to display signif- 
icant anomalies. For example, it is common to investi- 
gate the physics of disordered media at a marginal state 
of connectivity in terms of the geometry of the spanning 
cluster at the percolation threshold [[| 0]. First, it is 
clearly an advantage to use the percolation model be- 
cause a comprehensive set of exactly- and numerically- 
calculated critical exponents is available to describe not 
only its geometrical features, but also its dynamical 
(transport) properties. Second, the application of this ge- 
ometrical paradigm can be consistently justified through 
"the critical path method" S , a powerful approximation 
that has been successfully used M to estimate transport 



properties (e.g., permeability and electrical conductivity) 
of disordered porous materials. Accordingly, the trans- 
port in disordered media with a broad distribution of 
conducting elements should be dominated by those re- 
gions where the conductances are larger than some crit- 
ical value. This value is the largest conductance such 
that the set of conductances above this threshold forms 
a network that preserves the global connectivity of the 
system. In percolation terminology, this is equivalent to 
the conducting spanning cluster. 

The extent to which the self-similar characteristic of 
the critical percolation geometry can modify the displace- 
ment process is unclear. Several studies have been de- 
voted to the investigation of the displacement process 
through percolation porous media at criticality |l(], p"l| . 
More recently JlJ, [13]], the dynamics of viscous dis- 
placement through percolation clusters has been inves- 
tigated in the limiting condition of unit viscosity ratio 
m = ' H2, where /ii and \xi are the viscosities of the in- 
jected and displaced fluids, respectively. In this situation, 
the displacement front can be modeled by tracer particles 
following the streamlines of the flow, and the correspond- 
ing distributions of shortest path and minimal traveling 
time closely obey a scaling Ansatz |l4|, Subsequently 
fL6|| , the dynamics of viscous penetration through two- 
dimensional critical percolation networks has been inves- 
tigated in the limiting case of a very large viscosity ratio, 
m — ► oo. The results from this study indicate that the 
distribution of breakthrough time follows the same scal- 
ing behavior observed for the case m = 1 |l2|, As 
a consequence, it has been suggested [jlfj) that the pro- 
cess of viscous displacement through critical percolation 
networks might constitute a single universality class, in- 
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dependent of m. 

The spatial distributions of porosity and permeability 
in real rocks are often close to random. However, one 
cannot assume that the nature of their morphological 
disorder is necessarily uncorrelated, i.e, the probability 
for a site to be occupied is independent of the occupancy 
of other sites. In fact, the permeability of some rock for- 
mations can be consistently high over extended regions 
of space and low over others, characterizing in this way 
a correlated structure ||. In the case of sandstone, for 
example, the permeability is not the result of an uncorre- 
lated random process. Sand deposition by moving water 
or wind (and other mechanisms of geological scale) nat- 
urally imposes its own kind of correlations. A suitable 
mathematical approach to represent the geometry of the 
pore spaces and the corresponding transport properties 
is correlated percolation H |(J This ap- 

proach has been successfully used to model permeability 
fluctuations and also to explain the scale dependence of 
hydrodynamic dispersion coefficients in real porous ma- 
terials j22j. 

Our aim here is to extend the investigation on the 
displacement dynamics between two fluids thr oug h two- 
dimensional percolation clusters at criticality jl2[ |l3| to 
the case where the pore space display long-range spatial 
correlations. We focus on the case of viscous penetration 
with two immiscible fluids of unit viscosity ratio (m = 1) 
to study the effect of long-range correlations on the dis- 
tributions of traveling length and minimal traveling time. 

The organization of the paper is as follows. In Sec- 
tion II, we present the mathematical model to simulate 
long-range spatial correlations and show some geometri- 
cal features of the pore structures generated by this tech- 
nique. We also describe the dynamical model to simulate 
the process of viscous displacement in porous media. We 
show the results in Section III and Section IV is discus- 
sion and summary. 



where O is the Heavyside function and the parameter (j> is 
chosen to produce a lattice at the percolation threshold. 
In Figs, la and lb we show typical backbones extracted 
from uncorrelated and correlated networks. Long-range 
correlations in site occupancy gives rise to variations in 
the structural characteristics of the conducting backbone 
p9| . To illustrate this effect quantitatively, Fig. 2 com- 
pares the fractal dimension of the conducting backbone 
calculated for uncorrelated and for correlated networks 
with 7 = 0.5. Indeed, the fractal dimension of the back- 
bone is significantly larger for the correlated case. 

For a given correlated network at criticality, we choose 
two sites A and B belonging to the infinite cluster and 
separated by a distance r. In oil recovery these repre- 
sent the injection and production wells (see Fig. 3). We 
then extract the percolation backbone between these two 
points. To model incompressible flow through this disor- 
dered system, we assume that the lattice sites have neg- 
ligible volume and the allocated bonds are homogeneous 
elementary units of a porous material with constant per- 
meability k and flow area a. We also consider that the 
dynamics of fluid displacement is governed by viscous 
forces and that m — 1 (the invading and displaced fluids 
have the same viscosity) . Under these conditions and due 
to the strictly convective nature of the penetration pro- 
cess, the velocity at each elementary unit can be modeled 
in terms of Darcy's law 



(3) 



where Pi and Pj are the values of pressure at sites i and 
j, respectively, and I is the length of the bond. Conser- 
vation at each site of the backbone leads to the following 
set of linear algebraic equations: 



ka 



(4) 



II. MODEL 

We start by describing the geometry of the disor- 
dered system studied here. Our basic model of a porous 
medium is a two-dimensional site percolation cluster at 
criticality H [7) modified to introduce correlations among 
the occupancy units |l7|, |TS , 19, 2(], pl[ . The correlations 
are induced by means of the Fourier filtering method 
(FFM), where a set of random variables u(r) is intro- 
duced following a power-law correlation function 



(u(r)u(r + R)) oc R~ 



[0 < 7 < 2], (1) 



where 7 = 2 is the uncorrelated case and 7 « cor- 
responds to the maximum correlation. The correlated 
variables u(r) are used to define the occupancy C( r ) 01 
the sites 



C(r) = 0[0-u(r)] 



(2) 



where qij is the volume flow rate through the bond and 
the summation is taken over all bonds connected to a 
node i that belongs to the cluster. As a macroscopic 
boundary condition, we impose a constant flow rate Q 
between the injecting point A and the extracting point 
B. In practice, we apply a unit pressure drop between 
wells A and B, and calculate the solution of Eq. (Q) in 
terms of the pressure field by means of a standard sub- 
routine for sparse matrices. Due to the linearity of the 
system, the computed velocities at each bond, Vij, can be 
rescaled to give a fixed total flow rate Q, independent of 
the distance between A and B, and the realization of the 
porous medium. This resembles more closely oil recov- 
ery processes where constant flow is maintained instead 
of constant pressure drop. 

To simulate the displacement of fluid through the per- 
colation backbone, we first note that, under the condi- 
tions of unit viscosity ratio (m = 1) and for a fixed pres- 
sure drop between the wells, the pressure field remains 
constant during the propagation of the invading front 
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through the percolation network. Another consequence 
of this simplifying assumption is that the front of invad- 
ing fluid in any bond (ij) of the lattice advances locally 
with a constant velocity equal to Vij . This situation can 
be expressed as 



dt 



dx 



= 



(5) 



where Fij(x,t) denotes the interface between invading 
and displaced fluids, t is time and x corresponds to the 
local longitudinal coordinate within each elementary unit 
(bond) of the porous material. Equation (||) expresses the 
fact that the physical system considered here is always 
and everywhere convective for any value of the imposed 
flow rate Q. This behavior is entirely analogous to the 
convective (non-diffusive) regime of hydrodynamical dis- 
persion (|, [2| , where the unsteady transport of a neutral 
tracer in a carrier fluid flowing through a porous mate- 
rial is totally dominated by convective effects. In the 
absence of diffusive effects, the tracer samples the disor- 
dered medium by following the velocity streamlines. In 
the general case of hydrodynamical dispersion, however, 
diffusion might play a significant role. If the pore space is 
sufficiently heterogeneous, local zones of small velocities 
can be found, even under conditions of high overall flow 
rates. As a consequence, the propagation of the tracer 
front in these regions may be diffusion-like if the char- 
acteristic time for convection, t c = £/v, is greater than 
the typical diffusion time, = £ 2 / ' D m , where D m is the 
molecular diffusion coefficient of the tracer in the carrier 
fluid. 

Applying the analogy between fluid displacement and 
the convective propagation of a tracer through a disor- 
dered porous material, we can now put forward a random 
walk picture for the front penetration of the invading 
fluid. Here we follow Refs. (T^, [l3| and use the particle- 
launching algorithm (PLA) , where the movement of a set 
of (tracer) particles is statistically dictated by the local 
velocity field. In the PLA, each particle starting from the 
injection point A can travel through the medium along a 
different path connected to the recovery point B, taking 
steps of length I and duration tij = £/vij (Fig. 4). The 
probability pij that a tracer particle at node i selects an 
outgoing bond ij (a bond where > 0) is proportional 
to the velocity of flow on that bond, cx Vjj/ Y] k Vjk, 
where the summation on k is over all outgoing bonds. 



The traveling length I is the number of bonds present in 
path C. Among the ensemble of all paths {C}, we select 
the path C* that has the minimal traveling time, t m i n , 



t m in{C*) = mint(C) 



(7) 



This quantity corresponds to the breakthrough time of the 
displacing fluid. For a given realization of the percolation 
network, we compute all the traveling lengths and the 
minimal traveling time corresponding to the trajectories 
of 10 000 tracer particles. For a fixed value of r, this 
operation is repeated for 10 000 network realizations of 
size L x L, where L = 512, so L » r. We carried out 
simulations for different values of r and find that there 
is always a well-defined region where the distributions of 
P(l) and P(t m i n ) follow the scaling form |l4j 



P(z) = A z 



(8) 



where z denotes I or t m i n: z* is the maximum of the prob- 
ability distribution, the normalization constant is given 
by A z ~ (z*)^ 1 and the scaling function has the form 

rani 



f(y) = exp(-a z y **) . 
The exponents <p z and al z are related by 
<f>, = V(4 - 1) • 



(9) 



(10) 



Note that the scaling function / decreases sharply when 
z is smaller than z*. The lower cutoff is due to the con- 
straint, / > r. 

In Figs. 4(a) and 5(a) we show log-log plots of the 
probability densities P(l) and P(t m i n ), respectively, for 
five different values of "well" separation, r = 4, 8, 16, 32, 
and 64. For each curve, we determine the characteristic 
size z* as the peak of the distribution and plot z* on 
a double logarithmic scale. As shown in Figs. 4(b) and 
5(b), the results of our simulations indicate that both I* 
and t* min , respectively, have power-law dependences on 
the distance r, z* <~ r dz . The linear fit to the data yields 
the exponents d z for each distribution, namely, 



di = 1-13 ± 0.02 (correlated) 



(11) 



and 



III. RESULTS 

We investigate the effect of spatial long-range correla- 
tions on the distributions of traveling length and minimal 
traveling time. The traveling time t of a path C is defined 
as the sum of the time steps Uj through each bond ij be- 
longing to a connected path between A and B 

t = J2 *iJ ■ ( 6 ) 

(ij)ec 



<k = 1.75 ±0.03 (correlated). (12) 

The same exponents reported in |l^, [fj| for the case of 
flow through uncorrelated percolation networks (7 = 2.0) 
at constant flux are 

di ~ 1.21 (uncorrelated) (13) 

and 

dt ~ 1.33 (uncorrelated). (14) 
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Once more, the differences in these exponents for the cor- 
related and uncorrelated cases can be explained in terms 
of the morphology of the conducting backbone. As 7 de- 
creases, the backbone becomes gradually more compact 
p9| . This distinctive feature of the correlated geometry 
tends to reduce the value of di and augment the value of 
dt as the strength of the long-range correlations increases 
(i.e., 7 decreases). In the limiting case of a homoge- 
neous system, the corresponding exponents are e?; = 1 
and di = 2 §,|l| || . 

Figures 4(c) and 5(c) show the data collapse obtained 
by rescaling I and t m i n by their characteristic size, I* and 
t* min . Both scaled distributions are consistent with the 
scaling form of Eq. (JsJ) . From the least-square fit to the 
data in the scaling regions, we obtain the exponents 



and 



gi = 2.35 ± 0.05 {correlated) 



g t = 1.89 ± 0.04 (correlated). 



(15) 



(16) 



For uncorrelated pore networks subjected to the condi- 
tion of constant flux, the exponents are [O, Il3| 



gi ~ 2.0 (uncorrelated) 



and 



.9/ 



2.0 



(uncorrelated) . 



(17) 



(18) 



The differences between these distribution exponents 
have their origins in the different levels of compactness 
between correlated and uncorrelated clusters. These 
results are compatible with those of previous studies 
J[9[ p0| , which indicate that spatial correlations can 
change other critical exponents. 



IV. DISCUSSION 

The need for a better description of the geometrical 
features of the pore space has been the main conclu- 
sion of several recent experimental and theoretical studies 
on transport phenomena in disordered porous materials 



Q. It is therefore necessary to examine local aspects of 
the pore space morphology and relate them to the rele- 
vant mechanisms of momentum, heat and mass transfer 
in order to understand the important interplay between 
porous structure and phenomenology. From a concep- 
tual point of view, this task has been accomplished in 
many works, where computational simulations based on 
a detailed description of the pore space have been fairly 
successful in predicting and validating known correla- 
tions among transport properties of real porous media 

M I2I 
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|29fl . In the present work, we have in- 
vestigated the dynamics of immiscible fluid displacement 
using the framework of a percolation model for porous 
media that has been specially modified to introduce spa- 
tial long-range correlations among the occupancy units 
of permeability. This model is certainly a more real- 
istic description for the geometry of porous rocks and 
should lead to a better mathematical representation of 
their transport properties. 

Our results on the distributions of traveling length and 
minimal traveling time through correlated percolation 
networks show that spatial fluctuations in rock perme- 
ability can have significant consequences on the dynamics 
of fluid displacement. More precisely, we observed that 
the presence of long-range correlations can substantially 
modify the scaling exponents of these distributions and, 
therefore, their universality class. As in previous studies 
on the subject |2(J, we explain this change of behavior 
in terms of the morphological differences among uncorre- 
lated and correlated pore spaces generated at criticality. 
Compared to the uncorrelated structures, the backbone 
clusters of the correlated cluster has a more compact ge- 
ometry. The level of compactness depends, of course, on 
the degree of correlations introduced during the genera- 
tion process. Moreover, our results are consistent with 
the fact that the dynamical scaling exponents di and dt 
obtained for correlated geometries assume values inter- 
mediate between the uncorrelated and the homogeneous 
limiting cases. 
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FIG. 1: (a) A typical percolation lattice L — 256 for the 
uncorrelated case. Heavy lines correspond to the backbone, 
gray lines to dangling ends, and light gray lines to isolated 
"islands" (non-spanning clusters), (b) The same as in (a), 
but for the correlated case (7 = 0.5). 
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FIG. 2: Log-log plot of the backbone mass Mbb versus the grid 
size L for uncorrelated networks (squares) and correlated net- 
works (circles). The long-range correlated percolation struc- 
tures have been generated with 7 = 0.5. The solid lines are 
the least-square fits to the data with slopes corresponding to 
the fractal dimensions of the respective backbones, d bb . 




FIG. 3: The traveling paths of 10 000 tracers (L = 64, r = 16 
and 7 = 0.5). Heavy lines correspond to the bonds that re- 
ceive more than 6000 tracers, medium lines to those that re- 
ceive between 1 and 6000, and thin lines to those that receive 
none. 
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FIG. 4: (a) Log-log plot of traveling distance distribution 
P(l) for 7 = 0.5 and r = 4, 8, 16, 32, 64. (b) Log-log plot 
of the most probable value I* for traveling length versus the 
distance r. The straight line is the least-squares fit to the 
data, di — 1.13±0.02. (c) Data collapse obtained by rcscaling 
I with its characteristic value I* ~ r . The least-square fit 
to the data in the scaling region gives gi = 2.35 ± 0.05. 




FIG. 5: (a) Log-log plot of the minimum traveling time dis- 
tribution P(tmin) for 7 = 0.5 and r = 4, 8, 16, 32, 64. (b) Log- 
log plot of the most probable values for the minimal traveling 
time tmin versus the distance r. The straight line is the least- 
square fit to the data, with the number indicating the slope, 
d t — 1.75 ±0.03. (c) Data collapse obtained by rescaling t m in 
with its characteristic time t* min ~ r 1 ' 75 . The least-square fit 
to the data in the scaling region gives g t = 1.89 ± 0.04. 



